home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / procssng / ccs / ccs-11tl.lha / lbl / hips / sources / vfft / vfft.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-01-20  |  8.3 KB  |  375 lines
  1. /*
  2. %    VFFT .C
  3. %    Virtual-Very Fast Fourier Transform
  4. %
  5. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  6.  
  7. This software is copyright (C) by the Lawrence Berkeley Laboratory.
  8. Permission is granted to reproduce this software for non-commercial
  9. purposes provided that this notice is left intact.
  10.  
  11. It is acknowledged that the U.S. Government has rights to this software
  12. under Contract DE-AC03-765F00098 between the U.S.  Department of Energy
  13. and the University of California.
  14.  
  15. This software is provided as a professional and academic contribution
  16. for joint exchange. Thus, it is experimental, and is provided ``as is'',
  17. with no warranties of any kind whatsoever, no support, no promise of
  18. updates, or printed documentation. By using this software, you
  19. acknowledge that the Lawrence Berkeley Laboratory and Regents of the
  20. University of California shall have no liability with respect to the
  21. infringement of other copyrights by any part of this software.
  22.  
  23. For further information about this notice, contact William Johnston,
  24. Bld. 50B, Rm. 2239, Lawrence Berkeley Laboratory, Berkeley, CA, 94720.
  25. (wejohnston@lbl.gov)
  26.  
  27. For further information about this software, contact:
  28.     Jin Guojun
  29.     Bld. 50B, Rm. 2275, Lawrence Berkeley Laboratory, Berkeley, CA, 94720.
  30.     g_jin@lbl.gov
  31.  
  32. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  33. %
  34. %    Every block has its own write_header() before its processing.
  35. %
  36. % AUTHOR:    Guojun Jin - LBL    5/10/91
  37. */
  38. char    usage[]="options:\n\
  39. vfft [-double] [-D2] [-VV] [<] image [> VFFT]\n\
  40. \n\
  41. -double        output double VFFT    \n\
  42. -D2        2-dimension VFFT    \n\
  43. -VV        real and imaginary are in separated planes    \n";
  44.  
  45.  
  46. #include "complex.h"
  47. #include "header.def"
  48. #include "imagedef.h"
  49.  
  50. U_IMAGE    uimg;
  51.  
  52. MType    dimen1len, dimen2size, fsize, vsize;
  53.  
  54. #define    row    uimg.height
  55. #define    cln    uimg.width
  56. #define    frm    uimg.frames
  57.  
  58.  
  59. main(argc, argv)
  60. int    argc;
  61. char*    argv[];
  62. {
  63. bool    dflag=0, dimens=0, vflag=0;
  64. int    i, j, f, hrows, hcols, logfrms, logrows, logcols;
  65.  
  66. format_init(&uimg, IMAGE_INIT_TYPE, HIPS, -1, *argv, "S20-1");
  67.  
  68. for (i=1;i<argc;i++){
  69.     if (argv[i][0]=='-'){
  70.     switch(argv[i][1]) {
  71.     case 'd':
  72.         dflag++;    break;
  73.     case 'D':
  74.         dimens++;    break;
  75.     case 'V':
  76.         vflag++;    break;
  77.     default:
  78. errout:        usage_n_options(usage, i, argv[i]);
  79.     }
  80.     }
  81.     else    if ((in_fp=freopen(argv[i], "r", stdin)) == NULL)
  82.     syserr("can't open %s as input", argv[i]);
  83. }
  84.  
  85. io_test(fileno(in_fp), goto    errout);
  86.  
  87. (*uimg.header_handle)(HEADER_READ, &uimg, 0, 0);
  88.  
  89. if (uimg.in_form==IFMT_DOUBLE || uimg.in_form==IFMT_DBLCOM)
  90.     dflag++;
  91. if (uimg.in_form > IFMT_DBLCOM || uimg.in_form==IFMT_ASCII)
  92.     syserr("undesired input image format");
  93.  
  94. if (frm < 2) {
  95.     vflag = 0;
  96.     dimens = 1;
  97. }
  98. if (vflag)    dimens=0;
  99. fsize = row * cln;
  100. hrows = row >> 1;
  101. hcols = cln >> 1;
  102. dimen1len = hcols + 1;
  103. vsize = row * dimen1len;
  104.  
  105. logfrms=logrows=logcols= -1;
  106. for (i=0,j=1;i<12;i++,j+=j) {
  107.     if (j==row)    logrows=i;
  108.     if (j==cln)    logcols=i;
  109.     if (j==frm)    logfrms=i;
  110. }
  111.  
  112. if (logfrms == -1 || logrows == -1 || logcols == -1){
  113. mesg("be patient for non power of 2 processing\n");
  114. #include "realvfft.cxx"
  115. exit(0);
  116. }
  117.  
  118. w_init(MAX(MAX(hrows, hcols), 1<<logfrms-1));
  119. if (dflag){
  120. double    *temp,
  121.     *re_i = (double*)nzalloc(fsize, sizeof(*re_i)<<1, "ibuf"),
  122.     *im_i = re_i + fsize,
  123.     *re_o, *im_o, *cvt;
  124.  
  125.     if (!dimens){
  126.         re_o = (double*)nzalloc(frm*vsize, sizeof(*re_o)<<1, "obuf"),
  127.         im_o = re_o + frm*vsize;
  128.         cvt = re_i;
  129.     }
  130.     else    cvt = (double*)nzalloc(fsize, sizeof(*cvt)<<1, "cvt");
  131.  
  132.     if (vflag)
  133.         uimg.o_form = IFMT_DVVFFT3D;
  134.     else    uimg.o_form = IFMT_DVFFT3D + dimens;
  135.     uimg.pxl_out=16;
  136.     (*uimg.header_handle)(HEADER_WRITE, &uimg, argc, argv, True);
  137.  
  138.     for (f=0; f<frm; f++){
  139.         if (uimg.in_form == IFMT_DOUBLE)
  140.             temp = re_i;
  141.         else    temp = im_i;
  142.         upread(temp, uimg.pxl_in, fsize, in_fp);
  143.         switch(uimg.in_form){
  144.         case IFMT_BYTE:    btod(im_i, re_i, fsize);
  145.             break;
  146.         case IFMT_SHORT:    stod(im_i, re_i, fsize);
  147.             break;
  148.         case IFMT_LONG:    itod(im_i, re_i, fsize);
  149.             break;
  150.         case IFMT_FLOAT:    ftod(im_i, re_i, fsize);
  151.         }
  152.         for (i=0; i<fsize; i++)
  153.             im_i[i] = 0.;
  154.  
  155.         dvfft_2d(re_i, im_i, logrows, logcols);
  156.  
  157.         if (dimens){
  158.         register double    *re_in = re_i,    *im_in = im_i;
  159.             temp = cvt;
  160.             for (i=0; i<row; i++){
  161.                 for (j=0; j<dimen1len; j++){
  162.                 *temp++ = *re_in++;
  163.                 *temp++ = *im_in++;
  164.                 }
  165.                 re_in += cln - dimen1len;
  166.                 im_in += cln - dimen1len;
  167.             }
  168.             fwrite(cvt, sizeof(*cvt)<<1, vsize, out_fp);
  169.         }
  170.         else    {
  171.         register double    *re_c = re_o + f*vsize,
  172.                 *im_c = im_o + f*vsize,
  173.                 *re_in = re_i, *im_in = im_i;
  174.             for (i=0; i<row; i++){
  175.                 for (j=0; j<dimen1len; j++){
  176.                 *re_c++ = *re_in++;
  177.                 *im_c++ = *im_in++;
  178.                 }
  179.                 re_in += hcols-1;
  180.                 im_in += hcols-1;
  181.             }
  182.         }
  183.     }
  184.     if (!dimens){
  185.         w_load(1<<logfrms);
  186.         for (i=0; i<vsize; i++)
  187.         dvfftn(re_o+i, im_o+i, logfrms, vsize);
  188.         if (vflag)
  189.         for (f=0; f<frm; f++){
  190.             fwrite(re_o+f*vsize, sizeof(*re_o), vsize, out_fp);
  191.             fwrite(im_o+f*vsize, sizeof(*im_o), vsize, out_fp);
  192.         }
  193.         else{
  194.         register double    *fop,
  195.             *re_c = re_o,
  196.             *im_c = im_o;
  197.         for (f=0; f<frm; f++){
  198.             for (i=vsize, fop=cvt; i--;){
  199.             *fop++ = *re_c++;
  200.             *fop++ = *im_c++;
  201.             }
  202.             i = fwrite(cvt, sizeof(*cvt)<<1, vsize, out_fp);
  203.             if (i != vsize)
  204.             syserr("f%d [%d] write %d", f, vsize, i);
  205.         }
  206.         }
  207.     }
  208. }
  209. else    {
  210. float    *temp,
  211.     *re_i = (float*)nzalloc(fsize, sizeof(*re_i)<<1, "ibuf"),
  212.     *im_i = re_i + fsize,
  213.     *re_o, *im_o, *cvt;
  214.  
  215.     if (!dimens){
  216.         re_o = (float*)nzalloc(frm*vsize, sizeof(*re_o)<<1, "obuf"),
  217.         im_o = re_o + frm*vsize;
  218.         cvt = re_i;
  219.     }
  220.     else    cvt = (float*)nzalloc(fsize, sizeof(*cvt)<<1, "cvt");
  221.  
  222.     if (vflag)
  223.         uimg.o_form = IFMT_VVFFT3D;
  224.     else    uimg.o_form = IFMT_VFFT3D + dimens;
  225.     uimg.pxl_out = 8;
  226.     (*uimg.header_handle)(HEADER_WRITE, &uimg, argc, argv, True);
  227.  
  228.     for (f=0; f<frm; f++){
  229.         if (uimg.in_form == IFMT_FLOAT)
  230.             temp = re_i;
  231.         else    temp = im_i;
  232.         upread(temp, uimg.pxl_in, fsize, in_fp);
  233.         switch(uimg.in_form){
  234.         case IFMT_BYTE:    btof(im_i, re_i, fsize);
  235.             break;
  236.         case IFMT_SHORT:    stof(im_i, re_i, fsize);
  237.             break;
  238.         case IFMT_LONG:    itof(im_i, re_i, fsize);
  239.             break;
  240.         }
  241.         for (i=0; i<fsize; i++)
  242.             im_i[i] = 0.;
  243.  
  244.         vfft_2d(re_i, im_i, logrows, logcols);
  245.  
  246.         if (dimens){
  247.         register float    *re_in = re_i,    *im_in = im_i;
  248.             temp = cvt;
  249.             for (i=0; i<row; i++){
  250.                 for (j=0; j<dimen1len; j++){
  251.                 *temp++ = *re_in++;
  252.                 *temp++ = *im_in++;
  253.                 }
  254.                 re_in += cln - dimen1len;
  255.                 im_in += cln - dimen1len;
  256.             }
  257.             fwrite(cvt, sizeof(*cvt)<<1, vsize, out_fp);
  258.         }
  259.         else    {
  260.         register float    *re_c = re_o + f*vsize,
  261.                 *im_c = im_o + f*vsize,
  262.                 *re_in = re_i, *im_in = im_i;
  263.             for (i=0; i<row; i++){
  264.                 for (j=0; j<dimen1len; j++){
  265.                 *re_c++ = *re_in++;
  266.                 *im_c++ = *im_in++;
  267.                 }
  268.                 re_in += hcols-1;
  269.                 im_in += hcols-1;
  270.             }
  271.         }
  272.     }
  273.     if (!dimens){
  274.         w_load(1<<logfrms);
  275.         for (i=0; i<vsize; i++)
  276.         vfftn(re_o+i, im_o+i, logfrms, vsize);
  277.         if (vflag)
  278.         for (f=0; f<frm; f++){
  279.             fwrite(re_o+f*vsize, sizeof(*re_o), vsize, out_fp);
  280.             fwrite(im_o+f*vsize, sizeof(*im_o), vsize, out_fp);
  281.         }
  282.         else{
  283.         register float    *fop,
  284.             *re_c = re_o,
  285.             *im_c = im_o;
  286.         for (f=0; f<frm; f++){
  287.             for (i=vsize, fop=cvt; i--;){
  288.             *fop++ = *re_c++;
  289.             *fop++ = *im_c++;
  290.             }
  291.             i = fwrite(cvt, sizeof(*cvt)<<1, vsize, out_fp);
  292.             if (i != vsize)
  293.             syserr("f%d [%d] write %d", f, vsize, i);
  294.         }
  295.         }
  296.     }
  297. }
  298. }
  299.  
  300. btof(ibp, obp, n)
  301. register byte*    ibp;
  302. register float    *obp;
  303. register int    n;
  304. {
  305. register int    i;
  306. for (i=0; i<n; i++)    *obp++ = *ibp++;
  307. }
  308.  
  309. stof(ibp, obp, n)
  310. register unsigned short    *ibp;
  311. register float    *obp;
  312. register int    n;
  313. {
  314. register int    i;
  315. #ifdef    _DEBUG_
  316. message("%d short to float\n", n);
  317. #endif
  318. for (i=0; i<n; i++)    *obp++ = *ibp++;
  319. }
  320.  
  321. itof(ibp, obp, n)
  322. register int*    ibp, n;
  323. register float*    obp;
  324. {
  325. register int    i;
  326. #ifdef    _DEBUG_
  327. message("%d itof\n", n);
  328. #endif
  329. for (i=0; i<n; i++)    *obp++ = *ibp++;
  330. }
  331.  
  332. btod(ibp, obp, n)
  333. register byte*    ibp;
  334. register double    *obp;
  335. register int    n;
  336. {
  337. register int    i;
  338. for (i=0; i<n; i++)    *obp++ = *ibp++;
  339. }
  340.  
  341. stod(ibp, obp, n)
  342. register unsigned short    *ibp;
  343. register double    *obp;
  344. register int    n;
  345. {
  346. register int    i;
  347. #ifdef    _DEBUG_
  348. message("%d short to double\n", n);
  349. #endif
  350. for (i=0; i<n; i++)    *obp++ = *ibp++;
  351. }
  352.  
  353. itod(ibp, obp, n)
  354. register int*    ibp, n;
  355. register double    *obp;
  356. {
  357. register int    i;
  358. #ifdef    _DEBUG_
  359. message("%d itod\n", n);
  360. #endif
  361. for (i=0; i<n; i++)    *obp++ = *ibp++;
  362. }
  363.  
  364. ftod(ibp, obp, n)
  365. register float*    ibp;
  366. register double    *obp;
  367. register int    n;
  368. {
  369. register int    i;
  370. #ifdef    _DEBUG_
  371. message("%d FtoD\n", n);
  372. #endif
  373. for (i=0; i<n; i++)    *obp++ = *ibp++;
  374. }
  375.